%% pde_MF
clc;
clear;

index = 1;

syms x y z nu nu_r C0 Ca Cd
switch index
    case 1
        u1 = y.^4+z.^2;
        u2 = z.^4+x.^2;
        u3 = x.^4+y.^2;
        w1 = sin(y)+z;
        w2 = sin(z)+x;
        w3 = sin(x)+y;
        p = (2*x-1).*(2*y-1).*(2*z-1);
end

u1dx = diff(u1,x);
u1dy = diff(u1,y);
u1dz = diff(u1,z);
u1dxdx = diff(u1dx,x);
u1dydy = diff(u1dy,y);
u1dzdz = diff(u1dz,z);
u2dx = diff(u2,x);
u2dy = diff(u2,y);
u2dz = diff(u2,z);
u2dxdx = diff(u2dx,x);
u2dydy = diff(u2dy,y);
u2dzdz = diff(u2dz,z);
u3dx = diff(u3,x);
u3dy = diff(u3,y);
u3dz = diff(u3,z);
u3dxdx = diff(u3dx,x);
u3dydy = diff(u3dy,y);
u3dzdz = diff(u3dz,z);
w1dx = diff(w1,x);
w1dy = diff(w1,y);
w1dz = diff(w1,z);
w1dxdx = diff(w1dx,x);
w1dxdy = diff(w1dx,y);
w1dxdz = diff(w1dx,z);
w1dydy = diff(w1dy,y);
w1dzdz = diff(w1dz,z);
w2dx = diff(w2,x);
w2dy = diff(w2,y);
w2dz = diff(w2,z);
w2dxdx = diff(w2dx,x);
w2dydx = diff(w2dy,x);
w2dydy = diff(w2dy,y);
w2dydz = diff(w2dy,z);
w2dzdz = diff(w2dz,z);
w3dx = diff(w3,x);
w3dy = diff(w3,y);
w3dz = diff(w3,z);
w3dxdx = diff(w3dx,x);
w3dydy = diff(w3dy,y);
w3dzdx = diff(w3dz,x);
w3dzdy = diff(w3dz,y);
w3dzdz = diff(w3dz,z);
pdx = diff(p,x);
pdy = diff(p,y);
pdz = diff(p,z);

f11 = -(nu+nu_r)*(u1dxdx+u1dydy+u1dzdz) + (u1*u1dx+u2*u1dy+u3*u1dz) + pdx - 2*nu_r*(w3dy-w2dz);
f12 = -(nu+nu_r)*(u2dxdx+u2dydy+u2dzdz) + (u1*u2dx+u2*u2dy+u3*u2dz) + pdy - 2*nu_r*(-w3dx+w1dz);
f13 = -(nu+nu_r)*(u3dxdx+u3dydy+u3dzdz) + (u1*u3dx+u2*u3dy+u3*u3dz) + pdz - 2*nu_r*(w2dx-w1dy);
f21 = -(Ca+Cd)*(w1dxdx+w1dydy+w1dzdz) + (u1*w1dx+u2*w1dy+u3*w1dz) - (C0+Ca-Cd)*(w1dxdx+w2dydx+w3dzdx) + 4*nu_r*w1 - 2*nu_r*(u3dy-u2dz);
f22 = -(Ca+Cd)*(w2dxdx+w2dydy+w2dzdz) + (u1*w2dx+u2*w2dy+u3*w2dz) - (C0+Ca-Cd)*(w1dxdy+w2dydy+w3dzdy) + 4*nu_r*w2 - 2*nu_r*(-u3dx+u1dz);
f23 = -(Ca+Cd)*(w3dxdx+w3dydy+w3dzdz) + (u1*w3dx+u2*w3dy+u3*w3dz) - (C0+Ca-Cd)*(w1dxdz+w2dydz+w3dzdz) + 4*nu_r*w3 - 2*nu_r*(u2dx-u1dy);

disp("===u1导数===");
matlabFunction(simplify(u1dx), "Vars", {x,y,z})
matlabFunction(simplify(u1dy), "Vars", {x,y,z})
matlabFunction(simplify(u1dz), "Vars", {x,y,z})
disp("===u2导数===");
matlabFunction(simplify(u2dx), "Vars", {x,y,z})
matlabFunction(simplify(u2dy), "Vars", {x,y,z})
matlabFunction(simplify(u2dz), "Vars", {x,y,z})
disp("===u3导数===");
matlabFunction(simplify(u3dx), "Vars", {x,y,z})
matlabFunction(simplify(u3dy), "Vars", {x,y,z})
matlabFunction(simplify(u3dz), "Vars", {x,y,z})
disp("===w1导数===");
matlabFunction(simplify(w1dx), "Vars", {x,y,z})
matlabFunction(simplify(w1dy), "Vars", {x,y,z})
matlabFunction(simplify(w1dz), "Vars", {x,y,z})
disp("===w2导数===");
matlabFunction(simplify(w2dx), "Vars", {x,y,z})
matlabFunction(simplify(w2dy), "Vars", {x,y,z})
matlabFunction(simplify(w2dz), "Vars", {x,y,z})
disp("===w3导数===");
matlabFunction(simplify(w3dx), "Vars", {x,y,z})
matlabFunction(simplify(w3dy), "Vars", {x,y,z})
matlabFunction(simplify(w3dz), "Vars", {x,y,z})
disp("===p导数===");
matlabFunction(simplify(pdx), "Vars", {x,y,z})
matlabFunction(simplify(pdy), "Vars", {x,y,z})
matlabFunction(simplify(pdz), "Vars", {x,y,z})
disp("===定常MF方程===");
matlabFunction(simplify(f11), "Vars", {x,y,z,nu,nu_r,C0,Ca,Cd})
matlabFunction(simplify(f12), "Vars", {x,y,z,nu,nu_r,C0,Ca,Cd})
matlabFunction(simplify(f13), "Vars", {x,y,z,nu,nu_r,C0,Ca,Cd})
matlabFunction(simplify(f21), "Vars", {x,y,z,nu,nu_r,C0,Ca,Cd})
matlabFunction(simplify(f22), "Vars", {x,y,z,nu,nu_r,C0,Ca,Cd})
matlabFunction(simplify(f23), "Vars", {x,y,z,nu,nu_r,C0,Ca,Cd})